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Dense mixture of granules and liquid often shows a sever shear thickening and is called a dilatant 
fluid. We construct a fluid dynamics model for the dilatant fluid by introducing a phenomenological 
state variable for a local state of dispersed particles. With simple assumptions for an equation of the 
state variable, we demonstrate that the model can describe basic features of the dilatant fluid such 
as the stress-shear rate curve that represents discontinuous severe shear thickening, hysteresis upon 
changing shear rate, instantaneous hardening upon external impact. Analysis of the model reveals 
that the shear thickening fluid shows an instability in a shear flow for some regime and exhibits the 
shear thickening oscillation, i.e. the oscillatory shear flow alternating between the thickened and the 
relaxed states. Results of numerical simulations are presented for one and two-dimensional systems. 

PACS numbers: 83.80.Hj,83.60.Rs, 83.10.Ff,83.60.Wc 



I. INTRODUCTION 

One of the most common materials of the dilatant fluid 
is a dense mixture of cornstarch and water, and it can 
be used to demonstrate a number of counter-intuitive 
behaviors that the shear thickening medium shows: sud- 
den solidification upon externally applied stress, quick 
re-fluidization after removal of the stress, formation of 
holes and protrusions under strong vibration[ll, Q, etc. 

These behaviors come from severe shear thickening and 
hysteresis, that dense colloid or dense mixture of granules 
and liquid often show. The shear viscosity increases al- 
most discontinuously by orders of magnitude at a certain 
critical shear rateQ, which makes the fluid almost rigid 
against the sudden application of stress. It is called a "di- 
latant fluid" by analogy with the behavior of a granular 
medium when a granular medium is densely packed 
in a bag that is flexible but non-stretchable, it cannot 
be deformed because the volume is constant. The gran- 
ular medium must dilate upon deformation due to the 
principle of dilatancy by Reynolds^]. 

There are several peculiar features in the shear thick- 
ening of the dilatant fluid: (i) the thickening is so severe 
and instantaneous that it might be used even to make a 
body armor to stop a bullet [1] , (ii) the relaxation after re- 
moval of the external stress occurs within a few seconds, 
that is quick but not as instantaneous as in the thickening 
process, (iii) the medium in the thickened state behaves 
like a rigid material allowing little elastic deformation as 
long as it is under stress, (iv) the viscosity shows hystere- 
sis upon changing the shear rateQ, (v) noisy fluctuations 
have been observed in the response to an external shear 
stress in the thickening regime 0, ll]- 

Despite of the apparent analogy between the behav- 
iors shown by these media, it is not clear if the shear 



thickening of the dilatant fluid has something to do with 
the property of dilatancy of granular media. Origi- 
nally, the shear thickening in colloid systems were re- 
garded as a result of the disorder transition of the layer 
and/or string structure developed in the low shear rate 
regime [9l4l2|. The dispersed particles align due the 
shear flow to give shear thinning in low shear regime, 
but the turbulent motion in high shear regime destroys 
this structure to give shear thickening. Such layer 
and/or string structures have been observed in numer- 
ical simulations [l3| and experiments 14 1 , and in some 
cases the shear thickening occurs when the structure is 
broken 0. However, there are some other cases where no 
significant structure change are observed upon discon- 
tinuous shear t hickening 1 1 5l4 1 8j . Hydrocluster formation 
has been prop osed as an alternative origin of the shear 
thickening jl6l [TqI. [20j . Due to hydrodynamic interaction 
among particles in the fiuid, there is a certain condi- 
tion that clusters of particles grow and they can give 
large viscosity. Such a cluster structure of particles has 
been first identified in numerical simulations [l9i| . then 
suggested by SANS experiments [ll]. More direct obser- 



vation has been made using fast confocal microscopy [21 
Jamming is another possibility under active debate in 
recent years in connection with the glass transition. In 
dense granular system, the iamming can cause the diver- 
gence of viscosity H, d, [2^]26j . In connection with the 
dilatancy. Brown and Jaeger studied the discontinuous 
shear thickening and obtained somewhat empirical con- 
stitutive relations [27}. 

There are no microscopic theories for the dilatant fluid 
yet in the sense that the shear thickening is derived from 
the elementary interactions among constituents of the 
medium, i.e. granules and fluid, but there are a cou- 
ple of semi-empirical theories: the soft-glassy rheology 



(SGR) model|28| and the schematic mode couphng the- 
ory (MCT)[2j^The SGR model is the model based on 
the stochastic dynamics with the activation energy that 
depends on the stress. This model is extended to de- 
scribe the shear thickening by introducing the stress de- 
pendent effective temperature. MCT, which gives reason- 
able description of the glass transition, has been extended 
schematically by introducing shear rate dependent in- 
tegral kernel. Both of the theories are semi-empirical 
and have been demonstrated to show the discontinuous 
shear thickening, but they have not been incorporated in 
the fluid dynamics to study its flowing behavior of the 
medium. 

Recently, the present authors constructed a fluid dy- 
namic model for the dilatant fluid by phenomenologically 
introducing an internal state variable, which determines 
the viscosity of the medium (30|. The state variable it- 
self is determined by the local stress [Slj. The purpose 
of this paper is to present detailed study on the flow- 
ing property of the medium represented by the model. 
We demonstrate that the model shows the discontinu- 
ous shear thickening transition and the hysteresis upon 
changing the shear rate as has been observed in exper- 
iments. It is also shown that the steady shear flow be- 
comes unstable for a certain parameter range against the 
shear thickening oscillation, where the medium alternates 
between the thickened and the relaxed states. 

The paper is organized as follow. The model is in- 
troduced in Sec. 2, and it is examined for a simple uni- 
form shear flow configuration in Sec. 3. Similar analysis 
is given for the gravitational slope flow and Poiseuille 
flow in Sec. 4. The response to an impact is simulated 
in Sec. 5. Effects of inhomogeneity is studied in Sec. 6 by 
two-dimensional simulations. Summary and discussions 
are given in Sec. 7. 



II. MODEL 

The model is based on the fluid dynamics with an in- 
ternal state variable that describes the local structure 
of particles dispersed in the liquid. The viscosity of the 
medium is determined by the internal state, which in 
turn changes in response to the local shear stress. We 
introduce each element of the model in the following. 

a. Fluid dynamics: The dynamics of the medium as 
a fluid is represented by the velocity field v{r), and is 
governed by the hydrodynamic equation. 
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where the Lagrange derivative is introduced: 
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The symbols p, P, and aij represent the density, the 
pressure, and the (i, j) component of the viscous stress 




FIG. 1: Schematic pictures for granular configurations: a re- 
laxed state(a) and a jammed state(b). 



tensor a, respectively. The last term in Eq.(IT|) represents 
the body force on the fluid due to the gravitational accel- 
eration gi. We employ Einstein's rule for the summation 
over repeated suffixes. 

We consider the incompressible fluid, thus the pressure 
P is determined by the incompressible condition 



V -vir) =0. 



(3) 



The viscous stress tensor is assumed to be expressed 
through the ordinary relation 



with the shear rate tensor 

dvi dv-j 
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Note that Eq.Q does not represent a linear viscosity 
because the viscosity rj is not constant but depends on 
the internal state variable cf) of the medium. 

b. Internal state of the medium: The dilatant fluid 
contains dispersed granular particles, which provides the 
system with an internal degree of freedom for a macro- 
scopic description. FiglT] shows a schematic illustration 
for a relaxed state(a) and that for a jammed state(b). 
The internal state may have a vector or even higher or- 
der symmetry in general, but in this work we study a 
simple case where the state is represented by a scalar 
field (f>{r). We assign = for the relaxed state and 
= 1 for the jammed state. 

For a given flow field v{r), we assume that there exists 
a stationary value 0*, toward which the state variable (jj 
changes as 
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with the time scale t. 

We may assume that r is constant in the case where the 
internal state changes due to the thermal fluctuation or 
some other mechanism independent of the flowing field. 
However, we adopt the variable time scale r that is in- 
versely proportional to the local shear rate F, 



rV- 



(7) 
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with a dimcnsionlcss constant r, because it is more natu- 
ral to suppose that the state change is driven by the flow 
deformation. Note that this form of r docs not introduce 
a new time scale to the system and makes it respond 
quite peculiarly to an external impact. 

The stationary value 0* is determined by the local flow 
and we assume that it is an increasing function of the 
local stress S. We employ a simple form 



(t)*is) = 4>M 



1 + {s/s^) 



(8) 



with the characteristic shear stress 5*0. The parameter 
<j)M represents the value of the state variable in the high 
stress limit and should depend upon the volume frac- 
tion of the granules and some other parameters of the 
medium. 

For the scalar values of the shear rate F and the shear 
stress S in Eqs.Q and (O, we adopt the definitions 



Tr[7 • 7] 



5* = 



Tr[o- • a], 



(9) 



which reduce to the ordinary shear rate and shear stress 
in the case of simple shear flow. 

c. Viscosity: The shear thickening property of the 
model comes from the (/)-dependence of the viscosity, for 
which we assume 



■q{(j)) = 770 exp 



A- 



1 - • 



(10) 



with the viscosity in the relaxed state ryo and a dimension- 
less parameter A. We have introduced the Vogel-Fulcher 
type strong divergence at the jamming point (/) = 1 in or- 
der to represent severe thickening observed in the dilatant 
fluid. In Eq. pH]) . the state variable (j) plays an analogous 
role with the inverse temperature in the glass transition. 
Note that the state variable (j) cannot he cj) > 1 even 
when (/)M > 1 in Eq.® if we employ Eq.© because the 
shear rate vanishes F \ as (/) 1 due to the diverging 
viscosity. 

d. Unit system: For numerical presentation, we em- 
ploy the unit system where 



770 = 5*0 = p = 1, 



(11) 



namely, the time, length, and mass are measured by the 
units 



To 
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To, mo = pio, 



(12) 



respectively. The rate I/tq gives the scale for the shear 
rate where thickening occurs, and the length scale io is 
the corresponding hydrodynamic length scale. For the 
cornstarch suspension of 41 wt%0], these parameters 
may be estimated as Sq « 50 Pa, rjQ « lOPa-s. and 
p « lO'^kg/m'^; which give tq ss 0.2 s and £0 ss 5 cm. 
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FIG. 2: Simple flow configurations and the coordinate system: 
(a) shear flow, (b) gravitational slope flow, (c) PoiseuiUe flow, 
and (d) impact by a bullet. 



III. SIMPLE SHEAR FLOW UNDER 
EXTERNAL SHEAR STRESS 

First, we will study behaviors of the dilatant fluid for 
a simple shear flow under an externally applied shear 
stress (Fig. [D^a)). The velocity field is assumed to be 
V = {u{z,t), 0,0) and the external stress imposes the 
boundary condition 



S{z,t) 



z=±h 



(13) 



where we have introduced the notation for the shear 
stress 



5(z,i) = 77(0)7(z,i) 
with the shear rate 

du{z, t) 



dz 



(14) 
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h is the half width of the flow and Sc is the applied stress 
at the boundarics(Fig. [H^a)). Then, Eqs.(|T|) and dH) 
become 



P 



du{z, t) 



d 



dt dz 
dd^{z,t) 
dt 



S{z,t), (16) 
\^{z,t)\(4>{z,t)~c^,{S{z,t))), (17) 



In the following, we flrst examine a steady flow solu- 
tion, then perform the stability analysis for the solution 
and the numerical simulation for these equations of mo- 
tion. 
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FIG. 3: (Color online) The stress-shear rate relation for the 
viscosity given by Ea. (|18p for various (f>M with A = 1. The 
inset shows the plots in the logarithmic scale. 



Steady flow solution 



The steady solution for Eqs. ((I3|) 
obtained as 

c 

= <?!)*(S'o), 7 



P7)) can be readily 



7*(5'c). (18) 



From these equations, we can obtain the relationship 
between the stress and the shear rate, which is plotted 
in Figini for various values of (pM with A = 1. In the 
logarithmic plots, the straight line with the slope 1 cor- 
respond to the linear stress-shear rate relation with a 
constant differential viscosity. One can see that there are 
two regimes: the low viscosity regime in the low shear 
stress and the high viscosity regime in the high shear 
stress. Between the two regimes, there is a branch where 
the shear rate decreases for increasing shear stress. The 
state in the middle branch can be unstable against in- 
finitesimal perturbation. 

From this stress-shear rate relation, we expect there 
should be hysteresis upon changing the shear rate; if the 
system starts from the low shear rate on the lower branch, 
the stress increases continuously, but before the system 
reaches the end of the lower branch, it should jump to the 
upper branch by discontinuous increase of the stress. If 
the system starts from the high shear rate on the upper 
branch and the shear rate decreases, the stress should 
jump to the lower branch before the system reaches the 
end of the upper branch. This sudden increase/decrease 
of stress corresponds to the discontinuous change of vis- 
cosity in the shear thickening. 



2. Linear stability of the steady flow 



dependence is only on the z coordinate: 

vir,t) = {%z+du{z,t), 0,0), (j){r,t) = (j),{S^)+5(j){z,t), 

(19) 

then the dynamics is analyzed using Egs. p^ and (|17p . 
Even within this restriction, we will see the steady shear 
flow in the middle branch may become unstable and the 
oscillatory flow arises. 

The linearized equations for the perturbation are now 
given by 



d 

r — 

at 



92 



52 



(20) 



d>(z,<) = %(^ci):7jjj{z,t) + {-i+^:iij,)6ci){z,i^'^] 



where the primes denote the derivative by its argument, 
and we have introduced the abbreviated notations. 



dri{(j>) 



0=0. (So) 



d0*(S'c) 
dSc 



(22) 

Then, the growth rate ftk of the perturbation for the 
Fourier component with the wave number k in the z di- 
rection is determined by 



-7*(/>t 77*, rilk - 7*(-l -f (f>W*i*) 



= 0. 



(23) 



This gives a positive real part of flk for the wave num- 
ber k that satisfies 



< fc2 < kj. 
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(24) 



in the case d'j^/dSc < 0, i.e. Sc is in the unstable branch 
of the shear stress-shear rate curve. Since the smallest 
possible wave number k for the perturbation is 7r/(2/i) 
and ?7*/p is the kinematic viscosity for the steady flow, 
we can interpret this result in the way that the steady 
shear flow in the unstable branch is unstable as long as 
the system width is larger than the momentum diffusion 
length due to the viscosity. 

For a given external shear stress Sc in the unstable 
branch, the flow becomes unstable for the system wider 
than 2hc = n/kc, where the growth rate flk has a finite 
imaginary part Uc given by 



(25) 



Note that the scales of fee and uJc are typically set by 1 /£o 
and 1/to although their actual values depend on r and 
the other system parameters, i.e. A and (pM. 




Now, we examine the linear stability of the steady 
shear flow given by Ea. ([T8)) . For a full analysis, an arbi- 
trary perturbation should be allowed, but here we exam- 
ine the linear stability against the restricted perturbation 
where the velocity is in the z-direction and the spatial 



3. Shear thickening oscillation in unstable shear flow 

The oscillatory behavior of the shear flow in the un- 
stable regime can be seen by numerically integrating 



Time 

FIG. 4; (Color online) Oscillation of the average shear rate 
u{h)/h in the shear flow for h — 1.3, 2, and 3 with A — cj>M ~ 
1, r = 0.1, and = 1.1. The (green) line, that overlaps the 
plot for h = 1.3, shows the plot for f{t) = C1+C2 e^'''^ sm{ojt+ 
6) with cj = 4.0, r = 5.8, 9 = 0.86, ci = 0.33, and C2 = 0.056. 



FIG. 5: (Color online) Time development of (j}{z) (left) and 
u{z) (right) in the shear flow oscillation for A = (pM — 1, 
r — 0.1, Sc ~ 1.1, and h = 2. Only the positive parts of the 
flow {z > 0) are presented. 



A. Gravitational slope flow 



Eqs.dlSl) and dTT]) with Eqs.iO and In FigH the 

average shear rates u{h)/h for (anti-)symmetric solutions 
are plotted as a function of time for various system width 
h with the constant shear stress 5*0 = 1.1 in the unstable 
regime for A = 0m — 1 and r = 0.1. The initial state is 
prepared as the steady solution for 5*0 = 1. For this 
set of parameters, kc = 1.18, which gives he — 1.33 and 
ujc = 3.91. 

For h = 1.3, which is smaller than /ic, the flow shows 
overdumped sinusoidal oscillation with the angular fre- 
quency 4.0, which is close to ojc- For larger /i, the oscilla- 
tion becomes self-sustained and non-linear; the gradual 
buildups of the flow speed are followed by sudden drops. 

This non-linear oscillation of shear thickening fluid is 
shown in more detail in Fig|51 where the time develop- 
ment of the 4'{z) and u{z) are plotted. Only the positive 
half of the solution is plotted for the (anti-)symmetric 
solution. In the plots, the oscillation starts from the al- 
most uniform shear flow in the high viscosity state with a 
larger value of 0. This flow cannot be completely uniform 
because the external shear stress Sc is in the unstable 
branch, thus the flow speed builds up gradually as the 
internal state </> relaxes to reduce the viscosity, but even- 
tually, (f) starts increasing when the shear stress becomes 
large enough. Then, larger value of causes higher vis- 
cosity, which decelerates the flow speed, but this causes 
even higher value of (j) because the inertia stress due to 
the deceleration is added on the top of the stress by the 
shear flow, which results in the sudden drop of the flow 
speed. 



IV. GRAVITATIONAL SLOPE FLOW AND 
POISEUILLE FLOW 



Similar analyses are performed for a gravitational slope 
flow and Poiseuille pipe flow. 



For the gravitational slope flow (Figl^^b)), Eqs.([T]) and 
(pl should be solved with the boundary conditions 



I ' 



(7 • n , , s = 0, 

\z=h{x,y) 



(26) 



where we have assumed that the bottom of the flow is 
located at 2: = and the flow depth at (x, y) is given by 
/i(x, y); the vector n represents the normal vector to the 
flow surface. The gravitational body force is given by 

g^{g sin 0, 0, -5 cos 6) = (g,, , 0, -g^_) (27) 

with the slope angle 0. 

For the flow field v = (^(z, t), 0, 0), Eqs.dB) and dH) 
become 



du{z,t) dS{z,t) 



dt 







dz 

dPjz) 
dz 



(28) 
(29) 



50(M) 
dt 



|7(z,t)|(0(z,t)-0,(5(z,t))) (30) 



with the shear stress and the boundary conditions 
are given by 



i(0) = 0, 



du[z, t) 



dz 



= 0. 



(31) 



Eq. 



can be solved immediately to give the pressure 

P{z)=g^{h^z) + P^ (32) 



with the atmospheric pressure Pq. 

The steady solution for Eqs. (f28| and ()30|) under the 
boundary condition pip is given by 



(z) = 0*(5g(z)), 7(z) = 



Sg{z) 



(33) 



ri[(l>*{Sg{z)) 

with the gravitational shear stress 

Sg{z)=g\\{h~z). (34) 
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Gravitational Slope Flow 
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FIG. 6: (Color online) Steady Gravitational flows: (a) the 
flow speed proflles as a function of z, (b) the surface flow speed 
vs g\\ , and (c) the time development of the surface speed. 



FIG. 7: (Color online) Poiseuille flow: (a) the flow speed 
proflles as a function of r, (b) the flow flux vs pressure gradient 
AP/L, and (c) the time development of the flux. 



From these, the flow speed u(z) and the flux per unit 
width $G can be calculated by 



u{z) = / 7(2') dz' , 
Jo 



u{z)dz. (35) 



In FigEl the flow speed profiles and the surface speeds 
given by Eg. ([55)1 plotted for several sets of parameters. 
The depth dependences of the flow speed are shown in 
FiglH^a) for some values of g\\ ; For small g\\ , the flow speed 
depends upon the depth parabolically as in a Newtonian 
fluid, while, for larger the flow speed profile develops a 
convex part, which corresponds with the unstable branch 
of FiglS] in the shear flow. In FiglBJb) , the surface speed 
u{h) are plotted as a function of for some values of h. 
One can see they decreases for large , which means that 
the fluid flows slower for larger inclination angle. This is 
because the viscosity of the fluid becomes large in the 
high shear stress caused by large (/y . 



B. Poiseuille Flow 

Pressure driven pipe flow with the cylindrical symme- 
try around the x-axis (Figl2][c)) is governed by the equa- 
tion 



P 



du{r, t) 
dt 

dt 



(36) 



= -|7(r,t)|(0(r,t)-0,(5(r,t))) (37) 



with the shear stress and the shear rate 

du{r, t) 



S{r,t) =ri{(l))j{r,t), j{r,t) 



dr 



(38) 



Here, AP (> 0) is the pressure drop along the pipe over 
the length L, and r is the distance from the central axis: 
r = y^y^ + z^. 

The steady flow solution for this configuration is given 

by 



(r) = 045p(r)), 7(r) 



with the Poiseuille shear stress 

Sp{r) = --—r. 



Sp{r) 



ri(MSpir)) 



(39) 



(40) 



In Figl?! the flow speed proflles u{r) and the flow flux 
<!> defined as 



u{r)2TTr dr 



(41) 



are plotted. General features of the flow is analogous to 
those of the gravitational flow, and the flow flux decreases 
upon increasing the pressure gradient for the large pres- 
sure gradient because of the shear thickening. 



Shear thickening oscillation in gravitational flow 
and Poiseuille flow 



These steady flows become unstable when the shear 
stress is in the range of the unstable branch at some re- 
gion of the flow. The oscillations in the surface flow speed 
and the flow flux are plotted for the gravitational and 
Poiseuille flow in FigslHKc) and[7Jc), respectively The 
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shear thickening oscihation appears in a large enough 
system for a certain range of external drive or AP/L; 
The system length scale should be larger than the viscous 
length scale of the flow, and the external drive should be 
in the range where some part of the flow is in the unsta- 
ble branch. From the plots, one can see the oscillation 
disappears when the external drive is cither too small or 
too large. In the former case, the fluid behaves as New- 
tonian while in the latter case the size of the unstable 
region becomes too small. The shape of the oscillation 
in the non-linear oscillation regime is saw-teeth like, i.e. 
gradual increases followed by sudden drops, as we have 
discussed in the simple shear flow case. 

The spatial variation of oscillatory flow are shown in 
FiglH The general feature is the same with that of the 
shear flow, but the value of is zero at the surface of 
the gravitational flow and at the center of Poiscuillc flow 
because the shear is zero. 



V. RESPONSE TO AN EXTERNAL IMPACT 

One of the peculiar features of the dilatant fluid is in- 
stantaneous hardening by an external impact. It hardens 
almost immediately upon application of an external im- 
pact and allows little deformation like rigid material. It 
has been demonstrated that the hardening is so rapid 
that the material can be used for a body armor to stop 
a bullet Such instantaneous hardening cannot be ex- 
plained by the transformation between steady conflgu- 
rations of granules, but must be a result of the failure 
to rearrange the granular conflguration due to some ob- 
struction. Upon sudden impact, the granules are inhib- 
ited to rearrange their conflgurations due to to either dis- 
sipation by the interstitial fluid or the jamming by direct 
contacts. In the case of slow deformation, the stress is 
low and the lubrication due to the fluid allows the gran- 
ules to re-arrange themselves so that they can pass each 
other. 

In the present model, this aspect of the medium is rep- 
resented by Eq.(I7]) that the relaxation rate of the inter- 
nal state is proportional to the shear rate. For a sudden 
deformation, the state variable changes to a high stress 
value as the medium deforms; When the medium is dense 
{(pM ^ 1) and the external impact is strong enough, the 
state reaches the jammed state after a certain amount of 
deformation, which is almost independent of the speed of 
deformation. 

In order to demonstrate this aspect of the model, we 
perform simple simulations that the layer of fluid of the 
thickness h is driven by a sudden motion of the upper 
boundary wall at z = ft, with the flxed lower boundary at 
z = (FigUJi). Let U{t) = u{h,t) be the velocity of the 
upper wall. Initially, the fluid is at rest, 

u{z, t) = 0, (j){z, t) = 0, U{t) = for < < 0, (42) 

then the upper wall is moved suddenly by the velocity uq 
at t = 0, J7(0) = uq. For i > 0, the velocity of the upper 



wall is determined by 

dU{t) 

m ; — = - 

dt 



du{z, t) 



dz 



z—h 



(43) 



with Eas. (jl6p and (|17p . where m is the mass of the upper 
wall per unit length. 

Fig shows the displacement X{t) of the upper wall. 



X{t) = f U{t')dt', 
Jo 



(44) 



for the three cases, (pM ~ 0.8, 1, and 2 for various ini- 
tial speeds uq increases. The wall decelerates rapidly as 
the fluid thickens in response to the stress, and eventu- 
ally stops. For 0M = 0.8, the final wall displacement 
increases as the initial speed uq. On the other hand, for 
(pM — 1 and 2, the final displacement hardly depend on 
uq when uq > 5. This is because the fluid gets jammed 
at a certain strain as it deforms, and cannot deforms fur- 
ther. However, when the initial speed is small enough, 
the upper wall docs not stop quickly because the fluid 
does not thicken as shown FigElJc). 



VI. TWO DIMENSIONAL INHOMOGENEOUS 
FLOW 

Now, we present the results of numerical simulations 
for two dimensional system in the simple shear conflgu- 
ration (Fig|5Ja)) in order to examine how the inhomo- 
geneity in the x direction affects the system behavior, 
especially in the case of shear thickening oscillation. 

The velocity field is assumed to be in the x — z plane, 
V = {u{x, z,t)^0,w{x, z,t)), and in the x direction we 
employ the periodic boundary condition with the system 
length L. We take L = lOh in the present simulations. 

The fluid dynamic equation ([l} is integrated using the 
standard MAC(Marker-and-Cen) method[3i| for the in- 
compressible fluid, and \'V-v\ is kept less than 10~^°. Eu- 
ler method is employed for the time integration of Eqs.([T]) 
and ©. 

The motion of the plates at z ~ ±/i is controlled so 
that the average shear stress on the plate is equal to ^e. 



L 



dx = Sp 



z=±h 



(45) 



As for the initial configuration at i = 0, we assume 
that the fiuid is at rest and the state variable (p is close 
to zero with small fluctuations introduced at every com- 
putational grid point r^. 



t;(r,0)=0, ,^(n;,0) = e., 



(46) 



where is a random variable uniformly distributed over 
[0, e) with a small parameter e. We take e = 10~^. Note 
that, for the case of e = 0, all quantities do not depend 
on X, thus the simulations reduce to the one dimensional 
case given in Sec. III. 
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(a) Gravitational Siope Fiow (b) Poiseuiiie Fiow 




FIG. 8: (Color online) Time development of ^(z) (left) and u{z) (right) in the oscillatory flow of the gravitational slope flow 
(a) and the Poiseuiiie flow (b). The parameters are A = (pM = 1 and r = 0.1 with = 0.8 and h = 2 for the gravitational flow 
and with AP/L = 1 and R = 3 for Poiseuiiie flow. 
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FIG. 9: (Color online) The time dependence of the displace- 
ment X after the impact for the system of (pM = 0.8 (a), 1 
(b), and 2 (c) with the initial speed uq = 40, 20, 10, and 5, 
and for the system of (fiM = 1.0 with uo = 5, 2, and 1 (d). 
The other parameters are h = 2, r = 0.1, and A — m — 1. 



FIG. 10: (Color online) Flow diagram of the shear flow for 
= 0.85. The internal state variable (p does not depend on 
X in the region colored with gray and blue. The inset is the 
same diagram obtained for (pM = 1. The other parameters 
are A = 1, r = 0.1, with L = lOh. 



A. Flow diagrams 

Fig [To] shows a flow diagram in the Sc — h plane for 
4>M = 0.85 with A = 1 and r = 0.1 (the inset for (j>M = 
1). The diagrams are determined by the simulations at 
the points with marks. 

In the steady shear flow (grey) and the oscillatory flow 
(purple) regions, the initial fluctuations do not grow, thus 
the flows remain homogeneous in the x direction and are 
the same with those in the corresponding one dimen- 
sional cases in Sec. Ill (FigfTTj): the dashed lines show the 
boundary for the two regimes in the one-dimensional case 
given by 

kciS.) = ^ (47) 

using the definition of fee in Eg.ip^ as a function of 5*0. 
In the low 5*0 side, one can see that this coincides with 
the corresponding boundary in the two-dimensional case 
between the steady shear flow and the oscillatory flow. 

The major difference between the one and the two 
dimensional cases is that these two homogeneous flow 
regimes are limited to the smaller Sc side. In the larger 



region, the initial fluctuations in the state variable 4> 
grows, thus the flow results in the inhomogeneous flow 
in the case of — 0.85(pink) or the jammed flow in 
the case of (j>nj — 1.0 of the inset. In the following, we 
examine the flows in these two regimes. 



B. Inhomogeneous oscillatory flow 

First, we examine the flow for (j>M = 0.85. In this 
case, the viscosity does not diverge and the medium keeps 
flowing. In FigH^l the time evolution of the upper plate 
velocity Up is plotted along with the case without initial 
fluctuations. The flow shows irregular oscillation with 
smaller amplitude compared with the noiseless case. 

The snapshots of (j) for a single cycle of oscillation in 
FigHS] reveal that the whole system is not thickened and 
the oscillation is governed by a few thickening bands. At 
the time when Up reaches its minimum (a), the thickening 
branch with the high value of (j) (the red region) is being 
extended along the direction of (1,1). As the system 
flows, this branch is stretched and Up gradually increases 
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3.0 r 



0M = 1-0, h = 4,Se; = 0.75 
(jjM = 0.85, h^i,Se = 0.95 




(a) t = 4( 



time 

FIG. 11: (Color online) Time development of the upper 
plate velocity Up in the oscillatory flow regime in the two- 
dimensional simulation. The initial fluctuations decay quickly 
and the flows show homogeneous oscillation as in the case of 
the one-dimensional system. 




time 



FIG. 12: (Color online) The time evolution of the upper plate 
velocity Up in the inhomogeneous oscillatory flow regime with 
(j)M = 0.85. The other parameters are h = 5, L = 50, 5*0 — 
1.5, A — 1, and r = 0.1. The initial fluctuation is given by 
e = 10^**. The uniform oscillation flow with e = (the dashed 



line) is shown for comparison. 



(b) , and eventually the branch breaks off and Up reaches 
maximum (c). Then, high shear rate makes the broken 
branches extend again to the other side to cause sudden 
deceleration (d). The thickening branches first appear 
both in (1,1) and (1,-1) direction, but the latter tend 
to disappear and transforms into the (1,1) direction in 
the course of time, and only the thickening branches in 
the (1,1) direction remain. 



C. Jamming caused by the instability 

For the case of (j)M = 1, the viscosity can diverges 
and the instability of the homogeneous flow causes the 
jamming to stop the fiow. 

FigslTiland[T51shows the simulation results for <j>M ~ 1- 




FIG. 13: (Color online) The snapshots of the state variable 
(j> taken during a cycle of oscillation presented in Fig |12l The 
arrows indicate flow velocity. 




FIG. 14: The spatial variation of viscosity 77 at jz = at 
several times in the jamming regime. The parameters are 
ft = 3, L = 30, Se = 1.5, 0M = 1 with ^ = 1 and r = 0.1. 



The time evolution of the viscosity distribution at z = 
is shovirn in Fig |14l Initially, the viscosity is rather uni- 
form virith some fluctuations, but peak structure appears 
soon around t ~ S with a certain characteristic length 
scale. Some of the peaks grow sharply and the thicken- 
ing regions strongly localize(i > 4), then the system is 
jammed and the flow stops within a period of oscillation 
of the homogeneous case(Fig lTC)) . 

FigHni shows the color maps of the pressure P (a) and 
the state variable </> (b) at t = A. The fluctuation of 
first stands out just below the moving plates, but higher (p 
regions form band structure, and extend along the prin- 
cipal axis of the shear deformation, (1,1) and (1,-1). 
Some of the bands reach the upper plate from the lower 
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(a) p^wm 'r^mp = -w 




FIG. 15: (Color online) (a) The spatial distribution of the 
pressure P and (b) the internal state variable in the system 
presented in Fig |14l at t = 4. 




£ = 

— £ = 10-* 



5 10 5 10 

time time 

FIG. 16: (Color online) The time evolution of (a) the up- 
per plate velocity Up, and (b) the maximum viscosity in the 
system presented in Fig |14l The dashed lines represent the 
oscillatory flow without fluctuations. 

plate, and they jam the system. 

In FiglWa). we present the velocity of upper plate Up 
as a function of time. The solid line shows the time evo- 
lution starting from the internal state with fluctuations, 
and the dotted line represents the case without fluctua- 
tions for comparison. The state variable suddenly loses 
homogeneity at t = 2.3, then thickening branches ap- 
pears, and the velocity Up drops to zero. The maximum 
viscosity is always found inside the thickening branch 
for t > 2.3, and its value sharply increases as plotted 
in FigfTBlb). We cannot simulate the system up to the 
time when the plates motion actually stops because the 
numerical time integration becomes difficult as the vis- 
cosity becomes large, since it requires smaller time step. 
In the present case, however, we expect the system is 
jammed because the decrease of Up and the increase of 
maximum viscosity arc rapid and monotonic. 



VII. SUMMARY AND DISCUSSIONS 

The shear thickening shown by a dense mixture of 
granules and fluid has some peculiar features: (i) instan- 
taneous hardening, (ii) fast relaxation to flowing state, 
(iii) rigid thickened state, (iv) hysteretic thickening tran- 
sition, (v) oscillatory flowing behavior. We constructed 



a fluid dynamics model by introducing a phenomenolog- 
ical state variable and showed that the model can de- 
scribe these features. Especially, we demonstrated that 
the shear thickening oscillation appears in various shear 
flow conflgurations. 

a. Comparison with visco- elasticity: A visco-elastic 
fluid such as a polymer melt shows analogous behavior to 
the dilatant fluid; It behaves like a solid in a short time 
scale and like a fluid in a long time scale. This may be 
compared with that of the dilatant fluid, i.e. instanta- 
neous hardening in response to an external impact and 
fluidization after relaxation of the applied stress. How- 
ever, there are some important differences; the visco- 
elastic medium changes its behavior according to the ob- 
servation time scale, and it also allows large elastic de- 
formation even in a short time solid like behavior. On 
the other hand, the dilatant fluid changes its behavior 
according to the stress, i.e. it stays hardened while it is 
under the stress and starts flowing within a few second 
after the removal of the applied stress. The dilatant fluid 
allow little deformation even under large stress. 

h. Shear stress thickening: In constructing the 
model, we assume the fluid is shear-stress thickening, i.e. 
the viscosity depends upon the state variable (j), and the 
steady value of the state variable 0* is determined by 
the local stress as in Eq.([5]). It is instructive to see what 
would happen if we assume as a function of the shear 
rate 0*(7). In this case, the viscosity is directly given 
as a function of the shear rate, ''7(^*(7)), thus we should 
not have a discontinuous thickening unless we assume a 
discontinuity either in 4>*('^) or in ri((j)). 

Experimentally, the most direct evidence for the shear- 
stress thickening should be obtained by the observation 
that the pipe flow flux is not monotonically increasing as 
a function of the applied pressure gradient. 

c. State variable: Although we introduced the state 
variable (j) phenomcnologically, we suppose that the vari- 
able represents a certain microscopic property of the 
medium, such as contact numbers between grains, associ- 
ated with the restrictions against local rearrangement of 
granular configuration. The variable could be a vector or 
a tensor, but we examined the scalar case for simplicity. 
It is notable that even the scalar state variable produces 
the anisotropic stress chain like structure in the system 
as we have seen in the two-dimensional simulations. 

d. Jamming and response to impact: A remarkable 
feature of the dilatant fluid is that hardening response 
is so instantaneous that the medium can be used for a 
body armor to stop a bullet 0. We believe that such an 
instantaneous severe hardening cannot be explained by 
a transformation between two steady states, i.e. from 
the fluid state under low stress to the rigid state under 
high stress. Instead, it must be a result that the rapid 
rearrangement in the granular configuration is inhibited. 
There are two possible mechanisms that inhibit the gran- 
ular rearrangement in the densely packed medium: the 
Reynolds dilatancy and the formation of stress chains. 
The Reynolds dilatancy can inhibit rearrangement by 
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the fluid friction because tlic rearrangement sliould in- 
duce tlie strong interstitial fluid flow through granules in 
order to compensate the local volume change caused by 
the dilatancy. The stress chains through direct contacts 
between granules can be formed by the application exter- 
nal stress and prevent granules from being rearranged. 

In the present model, such an aspect of hardening is 
represented by the fact that the relaxational time scale 
for the state change is not constant but proportional to 
the shear rate ([7|). Then, the state variable (j> reaches a 
steady value 0* in a time scale where the strain changes 
by the amount r. Consequently, the medium with (pM ^ 
1 can deform only up to a certain strain around r by a 
hard impact. 

e. Shear thickening oscillation: One of interesting 
results of the present model is the shear thickening os- 
cillation. The steady shear flow is unstable when the 
flow is in the unstable branch of the shear stress-shear 
rate curve and is wide enough compared to the diffusion 
length scale by the viscosity. In this case, the flow shows 
the oscillation between the thick state in the high shear 
regime and the thin state in the low shear regime. In 
two dimensions, uniform oscillation appears only in the 
smaller external shear stress region, but some oscillatory 
behavior remains even when inhomogeneity develops in 
the flow. 

The shear thickening oscillation in the homogeneous 
flow looks similar to the stick-slip motion in a frictional 
system. However, there are some important differences; 
the stick-slip motion starts with the sudden accelera- 
tion caused by the slip weakening resistance under the 
constant speed driving through a mechanism with finite 
rigidity, while the shear thickening oscillation starts with 
the sudden deceleration caused by the shear thickening 
transition under the constant force driving. 

Aradian and Gates have also studied the dynam- 
ics of shear thickening fluids and found oscillatory 
behaviorsfl^. They focused, however, on the regime 
where the structural relaxation time is much larger than 
the fluid dynamical time scale, which is appropriate for 
liquid crystal systems. Consequently, the time scale of 
the oscillation is of the order or larger than the struc- 
tural relaxation time scale, therefore the dynamics is 
completely dissipative. On the other hand, in the present 
case, the structural relaxation time of the state variable 
is set by the shear rate and always comparable with the 
fluid dynamical time scale, thus the period of the oscilla- 
tion is determined by the fluid dynamics. It is also clear 
that the stress from the inertia plays an important role 
in the thickening phase in the oscillation. 

/. Experiments: As for the hysteresis, Deegan ob- 
served the hysteresis loop of the viscosity under the os- 
cillatory stress for the cornstarch-fluid system and con- 
sidered it to be the mechanism for persistent holes [s^). 
His system is thinner and in the regime of mild shear 
thickening in comparison with the system studied by Fall 
et al@. In the present work, we try to model rather se- 
vere shear thickening in choosing the functional form of 
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FIG. 17: (Color online) Hysteresis loops in the shear stress 
vs the shear rate in response to oscillatory shear stress for 
various amplitudes. The oscillatory stress S{t) — Sc sin(ajt) 
is applied to the upper plate located sd, z = h with the lower 
plate fixed at 2: = 0. The system parameters are h = 0.5, 
(pM = 0.6, r = 2, and A = 1. 

the viscosity Eq. ([T0| , but the present model qualitatively 
reproduces main feature of his data, by adjusting the pa- 
rameters for (f>M and r fFig |17l compared with Fig. 5 of 
[3^). For this set of the parameters, the initial noise de- 
cays, thus the one and two dimensional simulations give 
the same results. 

Regarding the shear thickening oscillation, we could 
not find any literature on the experiment which shows 
clear oscillation. This may partly because the system 
needs to be large enough, typically wider than £q or a 
few centimeters, and the inhomogeneity may develops in 
three dimensional system, which can obscure the oscilla- 
tion. The noisy fluctuation often observed in rheometer 
experiments under constant stress may be explained by 
this mechanism0,Q. Nevertheless, one may easily notice 
the oscillation around 10 Hz simply by pouring the dense 
water-starch mixture out of a container. We believe that 
this oscillation should be explained by the shear thicken- 
ing oscillation. We are planning experiments that allow 
quantitative comparison with our results. 

Although in the different physical context, the clear 
oscillatory flows |35 [ a long with a discontinuous transi- 
tion and hysteresis [36l |37| have been observed in the liq- 
uid crystal system that shows the shear thinning due to 
the state dependent viscosity. Such behavior could be 
also described using the phenomenological model like the 
present one. 

Chaotic dynamics has been observed in dilute aqueous 
solutions of a surfactant in the experiment under the con- 
stant shear rate in the shear thickening regime[38l[39l|. It 
was interpreted as the stick-slip transition between the 
two states of the fluid structure, thus physical relevance 
to the present instability is not clear, but we also found 
the chaotic dynamics in the present model in the case 
of the constant relaxation time r with the large system 
width. 
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